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Abstract 

We propose a two-dimensional model for the organization of sta- 
bilized microtubules driven by molecular motors in an unconfined ge- 
ometry. In this model two kinds of dynamics are competing. The 
first one is purely diffusive, with an interaction between the rotational 
degrees of freedom, while the second one is a local drive, dependent 
on microtubule polarity. As a result, there is a configuration depen- 
dent driving field. Applying a molecular field approximation, we are 
able to derive continuum equations. A study on the solutions of these 
equations shows nonequilibrium inhomogeneous steady states in vari- 
ous regions of the parameter space. The presence and stability of such 
self-organized states are investigated in terms of entropy production. 
Numerical simulations confirm our analytic results. 

1 Introduction 

Using a theoretical physics approach, we introduce and analyze a stochastic 
model inspired by the self assembly process of the cytoskeleton. 

The cytoskeleton can be seen as the "infrastructure" of eukaryotic cells, 
providing for both (dynamically evolving) spatial structure and internal 
transport processes that are fundamental for the cell itself and its role in a 
living organism (see 0). We focus on questions concerning the statistical 
mechanics of this particular biological system, mainly its ability to assem- 
ble in a variety of symmetry-breaking phases, which are commonly believed 
to have nonequilibrium nature and are considered important for cell mor- 
phogenesis. The nonequilibrium forces that give rise to these phases are 
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usually ascribed to both the so called dynamic instability of microtubules 
(associated with a confining geometry, see for example Q) and the action 
of molecular motors. Dynamic instability is a nonequilibrium polimeriza- 
tion/depolimerization process that enables microtubules or actin filaments 
to exert forces on confining surfaces such as the cell membrane. Molecular 
motors are a much studied family of proteins that are able to generate active 
motion on cytoskeletal filaments. 

In vitro and numerical experiments, sometimes called self-organization 
assays (||, Q, 0, H|), show that ensembles of microtubules and active 
motors can self-organize under proper conditions, and this ability is strictly 
linked to motor activity. In these assays, motors are typically found in solu- 
ble complexes and operate when two or more filaments are present; dynamic 
instability doesn't seem to be determinant. 

In this paper we concentrate on the role of motor proteins and do not 
consider dynamic instability. In a previous article ([^) we have shown the- 
oretically that in motility assays, where motor proteins are adsorbed on 
surfaces and the "gliding" process of single microtubules is observed, rota- 
tional diffusion [U leads to purely diffusive dynamics. 

This fact brought us to investigate the breaking of rotational symmetry, 
in the form of pattern formation, as a many body effect. That is, we consider 
the excluded volume interaction of many filaments. This is different from 
cooperation of motor proteins that is typically investigated in the situation 
of muscle contraction; (see ]l0[ and [^5| for interesting recent examples of 
stochastic modeling in this context). 

Concisely, our problem can be stated in terms of existence of inhomoge- 
neous nonequilibrium steady states. From this viewpoint, our model is on 
the same level as many other models of non-biological systems, such as a 
fluid with convecting flow, or a nonequilibrium chemical reactor (this is true 
as long as we are not modeling the self-regulatory processes that take part 
in cell morphogenesis, see Q). 

What makes our choice distinct and particularly interesting for statistical 
mechanics is that the generalized force which keeps the system far from 
equilibrium is not, in general, a global field or some boundary condition, 
as is usually found in the literature on far from equilibrium systems, but a 
local release of energy associated with transport. 

The model, which is a two-dimensional lattice spin system, is presented 
in section ^. 

In order to discover the relevant variables of the real system we explored 
different ways to implement its microscopic dynamics [|ll|] . In this paper we 
limit the discussion to the simplest case we could find showing significant 
results. 

This point needs some clarification. Our working hypothesis is that the 
relevant features of the system are the competition between diffusion and 
motor drive, together with excluded volume interactions. The spirit of our 
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investigation is not to reproduce in detail the mechanical features of the mi- 
croscopic system but to analyze the behavior of averaged quantities. Many 
different dynamics can reproduce the right averages (we follow the concept 
of universality in statistical mechanics), so our job is to find the essential 
features that a microscopic dynamics needs to have to reproduce the desired 
behavior. Of course this leaves many possibilities open to choose a micro- 
scopic dynamics and one has to be guided by considerations on experimental 
models. These considerations are extremely useful to give an intuitive grasp 
on the interpretation of microscopic dynamics. Nevertheless, the rigorous 
justification of our microscopic model has to be sought in the results it gives 
for macroscopic, averaged observables. 

All of the dynamics we developed are based on the competition between 
a diffusive process (when motors are detached from filaments or not active) 
and a driven-diffusive one (when motors are active), this means that they 
fall in the class of competing dynamics models ( |l^| , fL5| ). 

In the system presented in this paper, the driven process is realized 
in a way that is resemblant to B. Schmittmann and R.K.P. Zia's driven 
diffusive systems ( [|l2[ , |l3| ) . This case is most easily interpreted in terms of 
many body motility assays, so we will stick to this interpretation throughout 
the paper. Unfortunately, we are not aware of any experiments of this 
kind focused on self-organization, so we dedicate section ^ to discuss the 
possibilities of employing this microscopic dynamics for the experiment in 
reference || and the results we obtained with a different choice for the driven 
process, which seemes more sound for this case. 

In section [3|, starting from microscopic dynamics, we use a mean field 
approximation to derive a set of four discrete evolution equations for the 
system. We discuss the problems that arise when one considers the contin- 
uum limit of these equations in order to obtain a more at hand system of 
differential equations. In particular, the result of this limit is related to the 
characteristic times of the two competing dynamics. 

Once we have obtained a set of differential equations on TZ 2 , we look for 
steady states in different regimes. In the two limiting cases of dominant diffu- 
sion and dominant drive we find, respectively, the usual homogeneous Gibbs 
states and blocked phases; whereas, when the two dynamics are mixed, a 
new class of entropy-producing, phase-separated states emerges. By linear 
stability analysis, we find in parameter space the instability region of the 
homogeneous states. 

Finally, in section || these results are compared with those obtained by 
numerical simulation. In order to tackle the problem of entropy production 
numerically, we apply some of the machinery related to the fluctuation the- 
orem by Gallavotti and Cohen (reviewed in JTH]), that has recently been 
established for stochastic dynamics ( JlTj , |l8| , [|19|1 ) . 
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2 Microscopic Dynamics 



We imagine that motors are adsorbed on a flat surface and can push the fila- 
ments as in motility assays. Our aim is to find a transition to self-organized, 
inhomogeneous states. Experiments of this kind, though in principle possi- 
ble, have not been tried to our knowledge. 

On the other hand, in the self-organization experiments described in || 
and |7|], microtubule-associated motors are linked in bundles by some other 
proteins, and these bundles are in solution. Links of these experiments to 
our model will be discussed in section ^. 

We discretize space, time, and the orientation of the filaments. So we 
define a 2-dimensional square lattice A and we imagine that each lattice 
site is either empty or occupied by the center of mass of a filament. This 
corresponds to associate to each lattice site x G A a spin a x which takes 
the null value when the site is empty. The other possible values of the 
spins are determined by the discretization on the orientational degree of 
freedom. With these assumptions, microtubules are treated as rigid rods, 
neglecting dynamic instability. The minimal choice is that spins take values 
in {+1, —1, +i, — i, 0}, corresponding to the four fundamental orientations of 
the filaments and to the empty site (see fig. ||). Our results will show that 
this choice is significant. 




Figure 1: Conventions adopted for lattice and occupation variables. Com- 
plex values of the spins are merely adopted to make our computational life 
easier. Microtubule directions are along the axes. Spins are 45o to lattice 
bonds so that interaction with all four nearest neighbors is symmetric. 

The thermodynamic equilibrium properties of the system are determined 
entirely by its hamiltonian H. We take 
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H = J^aX + K^atal (1) 

x,y x,y 
n.n. n.n. 

This is just the most generic form that H can take in our case, provided 
that the interaction is nearest neighbor. 

Intuitively, the first coupling constant, J < 0, stands for an interaction 
between the directions of the rods (insensitive to their orientations), and 
mimics excluded volume effects on directions. K, on the other hand, is 
sensitive to the presence or absence of the filaments. 

The choice of a nearest neighbor hamiltonian is connected to the problem 
of the length of the filaments. Reasoning geometrically, if a is the lattice 
spacing, a filament of length L should interact with ^ sites. On the other 
hand, we can do the statistical mechanics of this system after a preliminary 
coarse graining procedure. This renders the hamiltonian nearest neighbors, 
while the constants J and K keep memory of the original scale, becoming 
functions of the (mean) length L of the filaments (see for example f23[). 
The spin of each site takes the meaning of a cluster spin, so even though 
we speak of filaments as they were in 1 : 1 correspondence with our spins, 
this is merely a convention, because one spin may already be the result of 
mesoscopic averaging procedures. 

Prom expression (jl|) it is easy to see that the statics of this model is 
equivalent to that of the 2-D Blume-Emery-Griffiths model with only two 



relevant coupling constants (see [2C] for a complete mean- field analysis of the 
phase diagram). In fact, the effective levels reduce to three (er 2 = ±1,0), 
two of which are 2-fold degenerate. They correspond to the two possible 
directions (regardless of orientation) of the filaments and their absence. 

Since external fields are absent, our hamiltonian produces two kinds of 
spatially isotropic states, which are associated respectively with the well- 
known liquid/gas and isotropic/nematic transitions (|20|]). Notice that the 
phases that the system can exhibit are related with the assumption on the 
discretization of the rotational degrees of freedom. A finer discretization 
could bring to a larger number of phases. Anyway, the latter would be all 
be spatially isotropic, causing little change in the problem of the transition 
to inhomogeneous states. 

We now consider the time-dependent properties of the model. The evo- 
lution algorithm is made of two branches, of which the first is a diffusion 
and the second a driven diffusive process. Both processes take place in pres- 
ence of the hamiltonian H and include two kinds of elementary Montecarlo 
moves: 

• Nearest neighbor spin exchange (only effective when one of the two 
sites is empty) a x <-> a x i x' £ n.n.(x) (Kawasaki dynamics). 

• Local orientation a x ^ — > a' x ^ (Glauber dynamics). 
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In the diffusive steps, when the interaction with motors is not active, 
these two kinds of moves correspond respectively to translational and rota- 
tional diffusion of a microtubule. 

The probability of accepting a move is a modified Metropolis 

1 AH, 
A = -(1 + tanh— ) 

Let's now "turn on" the interaction with motor proteins, and see what 
happens. What we assume is that motor drive will act on the translational 
diffusion of the centers of mass of the filaments, trasforming it in a driven 
diffusion, and that rotational diffusion will be inhibited. 

In order to model motor activity, we modify A in the same way that 
is commonly used in the context of driven diffusive systems p^]. That is, 
AH(C,C) -> AH(C,C) ± E(C) when configuration C is favorable for 
motor pushing forward or pulling back the filament (see figure ||). H can be 
taken to be the same as in (|l]). E(C) is the driving field, and is proportional 
to the work done by the motors. 

Z 

A 
/ \ 




Figure 2: Motor action. AH is shifted by —E for Kawasaki exchanges with 
sites A and B, identified by the arrow head, and by +E for exchanges with 
C and D. 

This way of modeling molecular motor driven dynamics contains a num- 
ber of hypotheses. First of all, we suppose that the motors will push the 
microtubules every time they can. Secondly, the motors are spread with con- 
stant density and push one filament at a time (this last choice is intuitive 
for motility assays; in section |6| it will be discussed for the case of self- 
organization assays). Lastly, the motor action constrains the microtubules 
to preserve their orientations. The justification for this last hypothesis for 
a motility assay lies in the fact that as soon as two motors are attached to 
one filament, its direction is frozen except for the elastic fluctuations of its 
head and tail (see |26j and @). 

The fact that the pushing of the motor is actually a diagonal translation 
may cause some perplexities. On the other hand, our dynamics is set up 
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taking into account the fact that we have two degrees of freedom, transla- 
tional and rotational, that are independent. So they have to be activated in 
different elementary moves. The overall motion springs from a sequence of 
these moves. In this wiew, motor action is just a bias in the translational 
dynamics. We also tested (111]]) numerically a dynamics where motors push 
linearly (exchange is favoured with site Z in figure |2|) obtaining no quali- 
tative change in the results (but this choice complicates our analytic mean 
field calculations because of the fact that translations involve next nearest 
neighbours). 

We would like to make a final remark on the boundary conditions. We 
adopted periodic boundary conditions in our simulations. 

The simultaneous presence of two competing dynamics, each one with its 



own characteristic time, makes it reasonable to expect (see [15|) that possible 
steady states have a nonequilibrium (non Gibbsian) nature independently on 
the boundary conditions. This fact is confirmed by our mean field analysis. 

When considering the dynamics of active motors alone, the condition 
above ceases to be valid and we are in the case of a four species driven 
system. One should keep in mind the following remark. 

When £^0, the resulting driving field depends on the spatiotemporal 
configuration of the system. As it arises from a modification of the Metropo- 
lis acceptation probability, each translational step of the microscopic dy- 
namics automatically satisfies detailed balance with respect to the modified 
energy difference, 

W(C -> C) _ P eq (C) _ 

_ J3(AH±E(C)) 

W(C -> C) P eq (C) 

However, in general, detailed balance is not satisfied for periodic bound- 
ary conditions ]l2|. For this reason, one can expect to observe inhomo- 
geneous ("blocked") steady states even in the presence of the sole motor 
dynamics. 

To sum up, these are the main features of the model: 

It is a model with competing dynamics. Thus, the characteristic time 
scale ratio of the two dynamics plays an important role in its behavior |l4f| . 

It is related to so-called multiple species driven systems (113], p3| , |2l]), 
to which it reduces when Glauber-like dynamics (rotational diffusion) is not 
active and periodic boundary conditions are set. 

It can be interpreted as a reaction/diffusion model, in which four diffus- 
ing species, driven in four orthogonal directions, are subject to (equilibrium) 
chemical reactions controlled by the Ising hamiltonian ([!]). 
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3 Mean field approximation and continuum equa- 
tions 



The dynamics we described in the preceding section defines a Markov process 
on the lattice A. The time-dependent probability of a configuration C = 
Wx}x&A is given by the master equation (p^j, OX 



j t P(C, t) = (r p L p + T a L a )P(C, t) (2) 

where r p , r a are the characteristic times of the two dynamics. Physically, 
r a and rp represent a slow process with respect to the entire process of 
interaction between a motor and a filament. Therefore at the time scale we 
adopt, the dynamics is regarded as diffusive or driven respectively when a 
filament encounters on average a small or a large number of motors. 

Jj p , Jj a are the generators of the evolution when the motors are passive 
and active respectively. These operators contain the transition probabilities 
W(C,C) from configuration C to C", which can be easily written down 
following the description of the model given in the preceding section. For 
diffusive dynamics regulated by hamiltonian (|l|) (switched off motors) 

W{x\C^C) = 

x E,{(1 " $W>)§(1 - <k,o)(l - S a r ax ) [1(1 - Th^)] + 
+ (^ +SJ o)^,o^ +g ^[|(l-Th^)]} 



In the case of driven dynamics the corrections are that Glauber transi- 
tions (rotations) are inhibited, and AH(C, C) -> AH(C, C) ± E(C) when 
the configuration is such that the motors can do work (fig. ||). 



Instead of postulating, as is usually done [ 12|| , some mesoscopic equations 
on the basis of symmetries and physical considerations, we start from micro- 
scopic dynamics and develop a local mean field approximation. The main 
assertion of this approximation is that the probability of a configuration 
factorizes as 

PM = n p( a -) 

x£\ 

where p(cr x ) is the most general single site measure: 

i I 

(I runs on all the possible values of the spin). 

With these assumptions, the mean value of a function of n spins can be 
written 

\-T \@X\ 1 ■ ■ ■ 1 &X n 

))= E F(I 1 ,...,I n )[p Il (x 1 )...p In (x n )] 
h,...,I„ 



8 



We use this approximation to obtain the evolution equations of the first 
moments. Such discrete equations are sums of the mean values of some 
quantities computed in four different sites of the lattice, reaching next-next 
nearest neighbors. Namely, given a lattice versor a, one has to consider the 
cluster containing both the nearest neighbors of x and those of x + a. We 
get expressions that involve the parameters 

• H(x) = (1 — cr^) (density of holes). 

• M{x) = (cr 2 ,) (quadrupole moment; interpreted as the mean value of 
the net filament direction). 

• g(x) + i/(x) = (cr x ), identified with vector field D(x) = (g(x), f(x)) 
(orientation of the filaments) . 

Using Taylor expansions we obtain the following equations: 

H = - div[— V-ff + JHMVM - KH(1 - H)VH] + 

- rEdiv(HD) ^ ' 

M = div[(F) 2 V^ + JH(l - H)VM - KHMVH] + 

+ rEdw(HTD)] + (4) 
- (1 -r)(M + Th(4JM)) 

d\w[H 2 Vjj + JgHVM - KgHVH] + 
+ rEd x [(l - H + M)H] + (5) 

- (l-r) 5 (l + ±Th(4JM)) 

d\v[H 2 vjj - JfHVM - KfNVH] + 
- rEd y [(l -H-M)H] + (6) 

- (l-r)/(l-iTh(4JM)) 

Where T is the conjugate operator in R 2 . 

A continuum limit has been performed along an appropriate path in 
parameter space, so that one characteristic time and the lattice spacing 
a have been eliminated. Consequently, we are left with four significant 
(rescaled) parameters, J, K, E, r. The first two concern the equilibrium 
properties of the system, while the others are indicators of its being far from 
equilibrium. 

r gives the relative weight of the two dynamics and is related to the 
ratio r = ^ of the two characteristic times by the equality r = j^j. r is 
independent from the motor drive E; as matter of fact, E only appears in 
association with r in the equations, so that we interpret the quantity rE as 
the nonequilibrium generalized force generated by the motor action. 

In the equations above, r may be regarded as a free, phenomenological, 
parameter, as is not fixed by any assumption on our model. 
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Establishing the relationship between r and some underlying microscopic 
parameters, such as motor density and microtubule length, is a problem 
that falls beyond the reach of our model, and that requires a more detailed 
analysis of the physical system. 

For motility assays, the analysis of T.Duke and others (p6j) specifies 
the microscopic process as a Markovian stochastic evolution of the number 
n of motor proteins attached to a microtubule. Following this work, one 
can identify r a with the mean first passage time for going from n = 2 (two 
attached motors) to n = 1 (one attached motor), and r p with the mean time 
in which a microtubule is attached to or 1 motor, obtaining: 

L + 2d(P ( l L\ 
Ta= LTM^L\ ed ~- 1 -l) 



L 2d 
t p = - H 

V V 

where d is the mean distance covered by a filament before linking a 
motor; it can be related to the surface density of the motors and other 
physical parameters, v is the average speed of filaments, and L their length. 
This gives the expression for r 




We have already said (page ||) that microtubule length enters through the 
parameters of our model. The above formula is a way to relate the filament 
length to the parameter r via independent, more detailed modeling. 

Equations (||) ... © have the structure of reaction-diffusion equations, 
with no source for H, as filament number is fixed. The sources for M, g and 
/ arise from the interaction between filament directions. 

The existence of a source makes them different from other equations 
obtained in the same way for two-species driven diffusive systems ( |2l] . |27J ) • 
Another difference is that in this last case, having just two species and a 
field forcing in one spatial direction, mean field equations come out to be 
one dimensional, while ours are not. However, by the general consideration 
that patterns are usually lower-dimensional (see |2^]), we expect that our 
solutions will depend on one spatial variable, as well. In what follows we 
will give more specific justifications of this fact. 
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4 Steady states and stability 



A quick glance at our equations is sufficient to conclude that they always 
admit the homogeneous solution 

H = H 
f,9 = 

M = M which satisfies the Ising-like equation 
M + Th(4 JM) = 

where, by definition of the parameters, \M\ < H. 

This solution is the Gibbs equilibrium one when rE = 0. 

It is interesting to notice that, when E = 0, even if r ^ 0, our equations 
can be derived from an equilibrium free energy; that is, they have the form 
X = with F a suitable functional of the fields. In particular, equations 
(||) and (||) are decoupled from equations (||) and (Q). The first two describe 
a relaxation process to the equilibrium phases of our model, while the others 
admit as a unique stationary state the null field D = 0. These considerations 
cease to be true as soon as E ^ 0. 

In fact, when rE ^ 0, nonzero irreversible currents for / and g arise: 

J g = rEH{\ -H + M) (7) 

and similarly for Jr. These stationary currents are the mirror of a nonzero 
entropy production which we can measure numerically (see sec.[|), and, as 
with this last quantity, are linear in the nonequilibrium drive rE. 

We perform linear stability analysis around this homogeneous solution 
(resumed in figure ||), to find out that the system, far from equilibrium, for 
perturbations of short wavelength, becomes unstable in a region of the H, r 
plane. As can be seen from the picture, for low concentration of filaments 
and low r the system is always in a stable homogeneous state. One instability 
may arise after r reaches a critical value r c . 

The next thing to do is to verify that this instability leads to the inhomo- 
geneous solutions we expect. In doing this, we employ the same techniques 



used in [21] and [£7|], that is, successive substitutions leading to a nonlinear 
dynamical system for the variable H. Let's outline our procedure. By keep- 
ing terms to lowest order in the fields one has a hint at what the solutions 
look like. With this spirit, from equation (|3|) one obtains the relation 

VH = rEHD + ... (8) 

When r E (0,1) and the two dynamics are actually competing, Using (||) 
in eqn. (Q) we can write an expansion of M(x) around the Ising value M 

M = M + q[(1 - M)d 2 x H - (1 + M)dyH] + a 2 . . . (9) 
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eigenvalue 




- q 



q<q* 

1 




H 



H c 



Figure 3: Linear stability analysis results. Our methods coincide with those 
in Iffijj ]. Current terms are kept to lowest order, q is the wavelength of the 
perturbation. For_ q > q* the solution is always stable. For q < q* we find 
an instability if H 0-- 2H +^) -> I ~_ r _ j£ an( j ^ are f} xec ^ the stability 



2 I 



region looks like the one represented in the sketch above. 
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with 

1 

a = j 

( 1 - r )(l + SSE 4jBy) 

The last step is to use these two results in equations @ and @ . 

It is easy to acknowledge that the source terms for g and /, are equal to 
zero iff g and / themselves are. Thus, the two fields are tendentially brought 
to zero by the dynamics (equations (||) and (||)). 

However, the decay times of / and g are very different. Precisely, they 
are determined by the sign of M. What happens, then, is that one of the two, 
say /, goes rapidly to zero while the other doesn't, so equation (||) allows 
studying the stationary states that depend on just one spatial variable, and 
we can substitute relations (||) and (||) in equation (||) and obtain (similarly 
as in plj and ^?J) the ordinary differential equation for H 



/dH\2 



AH 4 



BH 3 + 2E 2 H 2 



r 2 E 2 (l + M)H- (1 -r)H\og\H\ (10) 



where A and B are constants of integration. The solutions of equation 
(10) belong to the class of cnoidal functions, and are investigated in the 
mathematical literature. These are periodic functions whose period and 
amplitude are related to each other, by some known expressions. If the 
system has a finite size, this constraint, toghether with minimization of 
energy, fixes the vavelength of an inhomogeneous state. 



H 



Figure 4: Plot of a solution of equation (fLG|) obtained with the boundary 
condition H — > 1 as x — > ±oo, which, given the total density, fixes the 
width of the stripe. 

Turning our attention to unbounded systems, we discover stripe-like 
states like that of figure |j. Our solutions are characterized by the fact 
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that H may reach the null value. So, if r is small enough, independently 
from the values of the coupling constants, there exist lines along which the 
density reaches its maximum, compatibly with excluded volume. Along the 
same lines we find a discontinuity jump of the orientation field D, whereas 
M fluctuates around the Ising value M modulated by the second deriva- 
tives of H (Notice that, at this stage, the constraint that M must always 
be lower than H has to be implemented; we verify that it is satisfied for 
suitable values of J). This behavior is confirmed by the results of our nu- 
merical simulations. For comparison, in figure || are shown some snapshots 
of stationary states taken from our simulations. 

Finally, we want to note that, for r = 1 (and E ^ 0), the four species are 
always simultaneously present and conserved, so the expansion around 
the Ising value of M ceases to be significant, and is substituted by the 
relation 

v#-**vi + ... 

In this formula the presence of the conjugate operator T imposes that 
the solutions depend on just one spatial variable, in agreement with a pure 
multi-species driven diffusive system, and we find so called "blocked" states 
present in the literature. 



5 Simulations. Entropy production 

We perform simulations at fixed total density equal to 1/2 on a 40 x 40 
square lattice. Montecarlo time unit r corresponds to 1600 moves, or one 
lattice sweep. E is typically one order of magnitude greater than J and K. 

Among the parameters we examine are the mean values of M and D, 
and the mean number of accepted rotational and translational Montecarlo 
moves (per sweep). These latter quantities measure the relative mobility of 
the system, that is, how much a state is "blocked" (see figures 

In order to calculate these means, data are sampled each 10 3 r. Total 
running time is about 10 6 r. The time to reach a steady state goes from 10 4 
to 10 5 t. Inhomogeneous steady states typically arise for big enough values 
of rE and have the form of a stripe either orthogonal or oblique with respect 
to the direction chosen by the filaments. An alternative state is shaped like 
a droplet (of arrows with the same direction). This last state is metastable 
for intermediate values of r, but its stability time diverges as r — > 1 (see 
figure |5[) . 

To analyze these states, we pay particular attention to the entropy pro- 
duction and the structure factor. 

Following Lebowitz and Spohn (||i~8|), we define the entropy production 
at the stationary state as 
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r = 



r = 0.7 



r = 0.8 



r = 0.8 



r= 1 



Figure 5: Steady states obtained with our simulations for J = 1, K = 1.5, 
E = 8 and different values of r. The case r = is the equilibrium state, 
and the system is not sensitive to the value of E. For r = 0.7 a stripe- 
like pattern is clearly identifiable. As r becomes closer to 1, droplet-like 
metastable states become more and more long lasting. For r = 1 we obtain 
blocked states. These results do not change sensibly as long as J > and 
J, K are significantly lower (about one order of magnitude) than E. K 
can be positive or negative. 
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Figure 6: Accepted translations per sweep at the steady state, mean value. 
The parameter is plot versus r, and E is fixed ( J = 1, K = 1.5, E = 8). 
The slope of the curve represents the mobility, which, after reaching the 
saturating value of zero has a jump to negative values, meaning that an 
increasing drive inhibits translations. The jump in this parameter corre- 
sponds to the critical value for r. Doubling E determines a change in the 
calculated points of less than four percent (this is also valid for figures [/], || 
and H) . 



p irr = p flux = ~(W(t)) 



with 



W{t) = E 



J bond,i(s)ds 

i=l bonds 



where Jbond,i is the net current of the i-th species along the bond in the 
driven dynamics, and the above integral is just a sum over Montecarlo times. 

With this definition, we discover two distinct regimes (fig. ||). 

In the first one, characterized by low r, a class of homogeneous entropy 
producing steady states is observable, and entropy production increases al- 
most linearly with r (see sec. ^|). 

In the second regime inhomogeneous steady states appear, and entropy 
production drops. The critical value of r that separates these two behaviors 
is also evident from the graphs of mean accepted moves per sweep (figures 
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Figure 7: Mean value of accepted rotations per sweep versus r at the steady 
state ( J = 1, K = 1.5, E = 8). A decided drop of this parameter is 
observable at critical r. 

^ and ^). We estimate this value to be r c ~ 0.45. Considering the dis- 
continuities in the parameters, we believe that the transition is of the first 
kind. 

On the other hand, the structure factor is defined as 

S(k) = (\G(k)\ 2 ) 

with 

L x 

In particular, we observe the values assumed by S(k) in correspondence 
with wave vectors k that are sensitive to stripe-like states in various direc- 
tions; namely k = (0, 1), (1,0), -^(1, 1), ^(1, —1)- If suitably normalized, 
these quantities take the value of unity when the stationary state is a stripe 
in a well-defined direction; typical values for droplet-like states are around 
1/2- 

Notice that time averaging of |G(fc)| 2 may be meaningless if the system 
explores a number of different inhomogeneous steady states, so one has to 
be very careful about the state being stable. 

Figure ^ is a plot of the structure factor of the steady state as a function 
of r. In the high r region, metastable states are subject to a slowing down of 
the dynamics and become long lasting. This is reflected by the two distinct 
behaviors of the structure factor that can be seen in the picture (see also 
figure H|). 
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Figure 8: Entropy production at the steady state as a function of r ( 
J = 1, K = 1.5, E = 8). After the parameter attains its maximum, there 
is a drop and the system begins to structure. In the high r region the 
situation is more confusing because of the appearance of metastable states. 
In fact droplets usually produce less entropy than stripes. When r = 1, 
entropy production of the steady states appears to be practically zero, even 
though the states are not absorbing, because the microscopic dynamics is, 
in principle, active. This effect is a consequence of the fact that the driving 
field E is very strong. 
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Figure 9: Structure factor as a function of r ( J = 1, isT = 1.5, E = 8) 
. This quantity starts to rise when r reaches its critical value. The dashed 
line represents metastable droplet-like states, which become long lasting as 
r increases, and stable for r = 1. The solid line stands for stripe-like states, 
independently on their orientations. 

To sum up, our numerical studies confirm the results of mean field anal- 
ysis, and reveal a more intricate phenomenology for higher values of r. Due 
to the appearance of a slowing down in the dynamics, we suspect the exis- 
tence of a second phase transition in this region, but more work is needed. 
In particular, a finite size scaling analysis could be important. 

6 Links to self-organization assays 

We proceed by suggesting a possible interpretation of the microscopic model 
in terms of self-organization assays (mainly ref. ||). 

In doing this there are a few main points to address. First of all, our 
spins are driven one at a time whereas filaments need to be at least two 
to be driven by the crosslinking soluble motor complexes of Q. Second, 
motor complexes do not have a fixed position, but can diffuse. Third, as the 
crosslinks between filaments can produce a torque, a new question arises: 
do we need to take into account a driven rotational dynamics ? 

Dealing with the first two points is quite a simple matter. 

In particular, the first affects the elementary move of the motors. We 
have tested a different dynamics for motor drive (|ll|) which is more intu- 
itively connectible with the experiments in question . Referring to figure ||, 
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the move of a motor is favoured by the field E(C) if the sites A and/or B 
are filled, and is an exchange with the next nearest neighbor Z. Numerical 
results show the emergence of inhomogeneous states and the same qualita- 
tive behavior as the ordinary model. New kinds of self-organized states, as 
the stripe parallel to filament direction of figure 10 appear. The continuum 
mean field theory becomes a bit more complicated, due to the appearance 
of terms in the powers of (1 — H) in our system of equations. These terms 
come from the need of a filled neighbor for a translation. Nevertheless, the 
solutions that we find remain valid in the limit of high filament density. 

As to the second point, at the scale we consider, motors act as a random 
field and it is not relevant whether this is a result of brownian motion or of 
a static configuration of the motors. 

The third question is more essential and involves the role of rotational 
driven dynamics in the process of self-organization. 

The fact that the forcing torque privilegiates the rotational moves align- 
ing the filaments does not affect the nature of the stationary states. In fact, 
in our ordinary model, inhomogeneous states arise thanks to the joint effect 
of the driven translations and the aligning due to the coupling J. So the 
torsional drive is not in competition with equilibrium thermodinamics and 
can be absorbed in the coupling constants J and K, provided the system ex- 
hibits orientational order. Only an analysis of the times of relaxation could 
distinguish between the two behaviors. 

Finally, the task of connecting the effective parameters of our model with 
the microscopic ones can be achieved with microscopic modeling of the kind 
of ref. [pfl], specifically oriented on self-organization assays. 



7 Conclusion 

We have introduced a model of nonequilibrium statistical mechanics whose 
most interesting feature is the fact that the generalized force is a dynamic, 
configuration dependent field. This feature has been inspired by the action 
of molecular motors on cytoskeletal filaments, so we have made quantitative 
effort throughout the paper to interpret our model in terms of systems that 
involve these objects. 

Through analytic (mean field) approach and simulations, we have found 
evidence for a nonequilibrium phase transition to inhomogeneous states. We 
were able to see that these inhomogeneous steady states are genuinely far 
from equilibrium, checking their microscopic currents and entropy produc- 
tion. This results may be seen as a theoretical evidence that a local drive 
in competition with diffusion and excluded volume effects are sufficient to 
reach self-organized states. 

Work in progress on some variations of the model described in this pa- 
per make us confident on the generality of this statement. We already 
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Figure 10: Stripe-like stationary state for the dynamics mimicking self- 
organization assays. The values of the parameters are J = 1, K = — 5, 
E = 2, r = 1/2. 

mentioned briefly two of these variations (pagg. [?] and |T^) . Other numerical 
experiments have been tried for quasi three dimensional geometries, that is, 
allowing filaments to cross each other when diffusing with translations, and 
different moves for motor action, all showing the emergence of inhomoge- 
neous stripe-like states. 

More work is needed to understand fully the statistical mechanical prop- 
erties of our model. In particular, we are working on three issues. 

First of all we want to perform a more careful analysis of the high r 
region, aiming to understand the nature of the stripe-like and droplet-like 
states and their relation with the finite size of the lattice. 

Secondly, we would like to know if the role played by the field M is the 
same when one has a less radical discretization of microtubule directions, 
for example if a £ S 1 . 

Lastly, the dynamical properties of the relaxation of this model to its 
stationary states are fully undiscovered. In particular, an analysis of this 
kind could be useful to understand the role played by torsional drive in 
dynamics resmbling self-organization assays. 
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